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Abstract 

The cluster-weighted model (CWM) is a mixture model with random covariates which allows 
for flexible clustering and density estimation of a random vector composed by a response variable 
and by a set of covariates. In this class of models, the generalized linear exponential CWM is here 
introduced especially for modeling bivariate data of mixed-type. Its natural counterpart, in the 
family of latent class models, is also deflned. Maximum likelihood parameter estimates are derived 
using the EM algorithm and model selection is carried out using the Bayesian information criterion 
(BIC). Artificial and real data are finally considered to exemplify and appreciate the proposed model. 

Key words; mixture models with random covariates; model-based clustering; cluster-weighted mod- 
els; generalized linear models. 



X 



1 Introduction 



Finite mixture models a re commonly employed in statistical modeling with two different purposes 



(JTitterington et al. 



19851 pp. 2-3). In indirect appl ications, they are used as semiparar rietric competitors 



of no npar amet ric density estimation tec hniques (see 



Titterington et al. 



200d p. 8 and 



Escobar and West 



19851 pp. 28-29, 



McLachlan and Peel 



19951) . On the other hand, in direct applications, finite mixture models 



are considered as a powerful device for clusteri ng by assuming that each mi xture component repre 



sents a group (or cluster) in the original data (see 



Fralev and Raftervl 



1998 and 



McLachlan and Basford 



1988). The areas o f application of mixture models is huge and range from biolo gy and medicine (see 



Schlattmann 



20091) to economics and marketing (see 



Wedel and Kamakura , 



20011) ; an overview is given 



in lMcLachlan and Peej |2000l) and iFriihwirth-Schnatteil pOOd ). 

This paper focuses both on direct and indirect apphcations. The framework is represented by data 
arising from a bivariate random vector (X,Y) , having joint density p{x,y), where Y is the response 
variable and X is the (random) covariate. A flexible fram e to analyze s uch data is represented by the 
family of mixt ure models with ran dom covariates (see, e.g., 



Hennig 



model (CWM; 



Gershenfeld 



20001 ) ■ to which the cluster- weighted 



1997|) belongs. The CWM principle consists in factorizing p{x,y), in each 
mixture component, into the product between the conditional density of Y\X = x and the marginal 
density of X by assuming a parametric funct fonal dependence of Y on x. Thus, CWMs embrace the 



potential of finite mixtures of regressions (see 



DeSarbo and Cronl. 



Chap ter 8) and of finite mixtures of distributions (see, 



1988 and 



Titterington et al. 



Friihwirth-Schnatter . 



1985 and 



2006, 



McLachlan and Peel 



20001 ): the idea of the former approach is adopted to model the conditional distribution of Y\x, while the 



principle of the latter is used to m odel both p (x, y) a. nd the marginal density of X 



The CWM was first proposed in 



Gershenfeldl (|l997^ . under the assumptions of local linear dependence 



and Gaussian component conditional/marginal densities (the so-called linear Gaussian CWM), in the 
context of media technology. Recently, the linear Gaussian CWM has been generalized from a twofold 
point of view: i) considering different component conditional/marginal distributi ons and ii) adopting dif - 



ferent functional relationships. As regards i), the CWM has been approached in 



Ingrassia et al. 



( 2012bh . 



in a statistical framework, including t distributions for the c omponent conditional/marginal densities. 



By using this model as a building block. 



Ingrassia et 



t CWMs for model-based clustering. As concern ii). 



al. 



(20 



Punzd ( 



3 ) pr esented a family of parsimonious linear 



20121 ) generalized the linear Gaussian CWM 



by allowing for possible nonlinear functional relationships as modeled by a polynomial curve. Finally, 
to allow the applic ability of the linear Gaussian CWM in high dimensional spaces of the covariates. 



Subedi et al 



(J2013f ) defined the (linear Gaussian) cluster-weighted factor analyzers which assumes a 
latent Gaussian factor structure for the covariates in each mixture component. 

So far, the applicability of the CWM seems to be restricte d to a real- valued rand om vector {X, Y) . A 



Ingrassia et al. 



( 2012a ) with the gen eralized 



McCuUagh and Nelder . 



19891 ) for the 



recent effort to solve this problem has been carried forward by[ 

linear Gaussian CWM which adopts a generalized linear model (see[ 

relationship of y on a; in each mixture component. This implies the possibility to model, for example, 

a count response via a Poisson distribution, a dichotomous response by a Bernoulli distribution, and a 

strictly positive response by a Gamma density. 

Starting from the generalized linear Gaussian CWM, the purpose of this paper is to make the CWM 
applicable for both direct and indirect applications on bivariate data of mixed type, that is data arising 
from two variables of different type (e.g., a count variable and a dichotomous variable). In the literature, 
little work has been done on clustering and distribution estimation of such data; recent papers are due to 



Hunt and JorgenserJ (J201lf ) and iBrowne and McNicholas J 2012 ). The standard m ixtu re-based approach 



is rep resented by the latent class model (LCM; see 



Vermunt and Magidson . 



2002 



and 



Hennig and Liao . 



20131 ): in our context, similarly to the CWM, the LCM fits p {x, y) by a mixture in which each component 
joint density is factorized as the product of the univariate densities of X and Y. This is made possible 
by assuming that X and Y are "locally" independent within each mixture component. However, this is 
a strong assumption in practical applications. On the contrary, our approach is based on extending the 
generalized linear Gaussian CWM to handle a possible no real- valued covariate; the extension is obtained 
by modeling the covariate, in each mixture component, with a distribution within the exponential family. 
The resulting model will be called generalized linear exponential CWM; it, differently from the standard 
LCM, allows for local linear dependence of y on X in each mixture component. 

The paper is organized as follows. In Section [2] the generalized linear exponential CWM is presented. 
Some notes about the relation with other mixture models are given in Section [3l In Section |4j maximum 
likelihood estimation of the model parameters is approached by considering, and detailing, the EM 
algorithm (see also Section [5]) . The method to select the number of components is given in Section [51 
Applications to artificial and real data are considered in Section [71 and the paper closes, with discussion 
and suggestions for further work, in Section [S] 



2 Genesis and formulation of the model 



2.1 The cluster- weighted model 



Let 



<{x,y;rp) = ^Trjp{x,y;^j 



(1) 



j=i 



be the finite mixture of distributions, with k components, used to estimate p{x,y). In ([T}, p {x,y;^j^ 
is the parametric distribution (with respect to ^,) associated to the jth component, ttj is the weight 
of the jth component, with ttj > and X]i'=i^j ~ l; ^'^'^ "^ = {'^'j^') j with n = (tti, . . . , Tr^) and 
^ = {^i, . . . ,^j,) , contains all the unknown parameters in the mixture. As usual, model ([ij implicitly 
assumes that the component joint distributions should all belong to the same parametric family. 
The mixture model ([T]) becomes a CWM when: 

• for each j, there is a parametric function 



E(Y\x,j;f3,) ^ fiYlj (x;f3 



3 y-^'f^jj 



(2) 



relating the expected value of the response Y to the covariate, where /3 ■ are regression parameters. 



In the following, for simplicity, we will consider ^ = (/?oj,/3ij) ; 
• the jth component joint distribution is factorized as 

P {x, y; ^j) = p [y\x; ^yu) p {^'^ ^x\j) , (3) 

where ^j — < ixij^irij p '^i^h ^y\j containing /S^. 
Summarizing, the CWM has equation 

k 

P (x, y;ip)^Y^ TTjp (^y\x; ^y\]) P (a;; ^x\j) ■ (4) 

2.2 Exponential family and "link" with generalized linear models 

If the distribution p (z) of a random variable Z is in the exponential family, then 

p(z;0,0)-exp|^^-^+c(z;0)|, (5) 

for specific functions a(-), b{-) and c(-). If is known, this is an exponential family with canonical 
parameter 9 which is function of the location parameter of the distribution. It may or may not be a 
two-parameter exponential family if is unknown. The function b{-) is of special importance because 
it describes the relationship between the mean value and the variance in the distribution (see, e.g.. 



McCullagh and Nelder . 



1989. pp. 28-29). In particular, E (Z) ^ fi ^ b' [9) and Var{Z) = 6" (0) a (0) 
where primes denote differentiation with respect to 9. The parameter (p is called the dispersion parameter 
and b" {9) is called the variance function. The variance function is often written as V (/i) = b" (0), where 
the notation V (/i) does not mean "the variance of /i" ; rather, V ( ii) indicates h ow the variance depends 
on /i in the distribution, and /i is in turn a function of 9 (see also 



Olsson . 



20021 pp. 37-40). 



It is well known that the exponential family model in ([5]) is strictly related with the generalized linear 
models. Here, a monotone and differentiable link function g (•) is introduced which relates the expected 
value fi, of the response Z, to the covariate X (a single covariate is enough for the scope of the present 
paper) through the relation 

gM = v = l3o + Pix. (6) 

Since the interest is now in the parameters (/3o, /3i) — (H, the distribution of Z given x will be denoted by 
p {z\x] f3, (f). The choice of the link function depends on the type of data. However, certain link functions 
are, in some sense, "natural" for certain distributions and they are called canonical links. In particular, 
the canonical link is the function g{-) such that <? (/^) — 9. Table [1] summarizes all the quantities seen 



so far for a few well-known distributions in the exponential family. It should be noted, however, that 
there is no guarantee that the canonical links will always provide the "best" model for a given set of 
data. In any particular application the data may exhibit peculiar behavior, or there may be theoretical 
justification for choosing links other than the canonical ones. For example, in the case of the Gamma 
distribution, the domain of the canonical link function is not the same as the permitted range of the 
mean. In particular, the linear predictor may be negative, which would give an impossible negative 
mean. When maximizing the likelihood, precautions must be taken to avoid this. An alternative is to 
use a noncanonical link function like the "log" . 

Table 1: Characteristics of some common distributions in the exponential family 





Gaussian 




Gamma 


Poisson 


Binomial 


Notation 


N{„,a^) 




g(m,^) 


P{i^) 


B(m,p) 


Support of Z 


(-00,00) 




(0,00) 


{0,1,...} 


{0,1,..., m} 


a(0) 


a^ 




i.-i 


1 


1 


6(.) 


02/2 




-\n{-e) 


exp(0) 


mln[l + exp(6')] 




1 

~2 

e 


i-+ln(2 
. 9 


K^) 


J/ In (vz) 
-0-1 


-ln(z!) 
exp (61) 


ln(T) 

exp (6*) 
"'l + exp(0) 


Canonical link 


identity 




^^-' 


log 


logit 


Var {Z) 


a2 




^,^/v 


M 


m/i(l-Ai) 


V{^^) 


^2 




^^ 


M 


mp(l-p) 



2.3 The generalized linear exponential CWM 

Consider the general formulation of a CWM in (|4]). By assuming that the distribution of X, given the 
component j, is within the exponential family ([5]), and by assuming a generalized linear model for Y on 
X for each j, we obtain the model 



p{x,y;ip) =^TTjp{y\x;f3j,(pY\j)p{x;0x\j,(f>x\j) , 

which will be referred to as the generalized linear exponential CWM. With respect to ^ we have ^ 
(/3j,0yb)' and ^^ I J = (0x|j, 0x|j)', whereby = U'y\i, ■ ■ ■ ,^Y\k) and |_,f = U'x\i, ■ ■ ■ ,^'} 



>X\k 



(7) 



Y\j - 

Note 



that the number d of free parameters in ([T]), that is the cardinality of ^p, depends on the exponential 
family distributions chosen as mixture components. 



3 Relationship with some common mixture-based approaches 



Let 



p{y\x;Tr,^Y) = ^■^jp{y\x;f3j,' 



V|jj , 



fc 



and 









(8) 



(9) 



(10) 



be, respectively, the (conditional) density of a finite mixture of generalized linear models of Y on x, the 
(marginal) density of a finite mixture of exponential family distributions for X, and the exponential- 
exponential LCM for [X, Y) where, in addition to the local independence assumption, an exponential 
family distribution is assumed for both X and Y in each group. In (jlOp. it results ^yi = {6y\j,4>y\j) 
and £y = I ^Y|i, • ■ • i^Fifc ) ■ To the author's knowledge, although the ratio underlying a LCM is well- 
known, its "exponential-exponential" variant has never been introduced. We have chosen to consider 
exponential family distributions in all the three models simply to relate them to our model. 

There is an important aspect which distinguishes the CWM in ([T]) from the mixture of generalized 
linear models in (|5]). The latter, in fact, belongs to the class of mi xture models with fixed covariate 
and, as such, does not consider modeling for X. As highlighted bv iHennieJ (|2000l ). this aspect makes 



the mixtures of regressions inadequate for most of the applications because, from a clustering point of 
view, they assume assignment independence: the probability for a point {x,y) to be generated by one 
of the k clusters depends on x only through Y , but the assignment of the data points to the clusters is 
independent of the covariate distribution. On the contrary, the CWM in ([7]) assumes a random covariate 
and allows for assignment dependence: the component covariate distributions can also be distinct (in 
terms of parameters) and they affect the assignment of the data points to the clusters. 

The next first two propositions give sufficient conditions such that the posterior probabilities related 
to models ([8]) and ([9]) coincide with those related to the generalized linear exponential CWM. Analogous 
re sults for the linear Gau s sian CWM, the l inear t CW M, and the po lynomial Gaussian CWM are given 



Ingrassia et al. 



I 2012bh 



Ingrassia et al. 



(|2013f ). and iPunzd (J2012n . respectively. The last proposition 



shows like model pop can be seen as nested in model ([7]). Proofs are given in the Appendix. 



Proposition 1. Given tv and ^y, if &x\i = ■ ■ ■ = Ox\k — ^x o,nd 4'x\i = ■ ■ ■ = 4'x\k 
([7]) generates the same posterior probabilities of component-membership of model ([5]). 



hx, then model 



Proposition 2. Given n and ^x: *//3i = ■ ■ ■ = /3fc = /9 o,nd (t>Y\i = ■ ■ ■ = (t>Y\k = 'Py, then model ([7]) 



generates the same posterior probabilities of component-membership of model ^ . 



Proposition 3. Model (|10p can be seen as a special case of model ([7]) when /3i 



I3ik = 0. 



We finally remark that, Proposition [T] highlights the cases where the generalized linear exponential 
CWM loses its property of assignment dependence in favor of assignment independence: the single 
covariate has the same (in terms of parameters) exponential family distribution in each group which, as 
such, can not affect the clustering results. On the contrary. Proposition [5] individuates a sort of "strong" 
assignment dependence for our model since clustering will only depend on X. Finally, Proposition [3] 
shows as model ([7]) generalize s its natural counterpart in the latent class modeling frame by relaxing the 



strong, and often unrealistic (jVermunt and Maeidson 



2002f) . assumption of local independence. 



4 The EM algorithm for maximum hkehhood estimation 

Given n observations (a;i, yi) , . . . , (xn, y„) from {X, Y) , the observed-data log-likelihood for the gen- 
eralized linear exponential CWM, when k is supposed to be pre-assigned, can be written as 



l{^)^Y.^n 



^TTjp{y,\xf, f3j ,(l)Y\j)p{xt;9x\j,(l>x\j) 



(11) 



The EM algorithm ([Dempster et all 119771) can be used to maximize I (ip) in order to find maximum 
likelihood (ML) estimates for the d unknown parameters of the generalized linear exponential CWM. 
Once k is assigned, it basically takes into account the complete-data log-likelihood 

n k n k n k 

i — l j — 1 2 — 1 J — 1 2—1 J — 1 

where z^ = 1 if {xi,yi) comes from component j and Zy = otherwise. The E and M steps of the 
algorithm can be detailed as follows. 



4.1 E-step 

The E-step, on the [q + l)th iteration, q = 0, 1, . . ., requires the calculation of the expectation of Ic (ip) 
given the observed data (a;i,yi) ,...,(x„,y„) and the provisional estimate ip , of ip, arising from 
the previous iteration. As Ic (xp) is linear in the unobservable data Zij, the E-step simply requires the 
calculation of the current conditional expectation of Zij given the observed sample, where Zij is the 



random variable corresponding to Zij. In particular, for i — 1, . . . ,n and j = 1, . . . , fc, it follows that 

i^,,(,) [Z^, |(x.,2/.)'] = -^+'^ = ^^ . '^\A '' '^^ (13) 



p[x^,y^;ip 



(9) 



which corresponds to the posterior probability that the unlabeled observation {xi, yi) belongs to the jth 
component of the mixture, using the current fit i/j^'^' for "0. 

4.2 M-step 

In the M-step, on the {q + l)th iteration, g = 0, 1, . . ., to maximize the expectation of Ic with respect to 
ip, the values z.y in (fT^ are simply replaced by their current expectations r|^ obtained in (J13p . This 
leads to 

n k n k n k 

i—l j — 1 i—l j — 1 2—1 J — 1 

(14) 
As the three terms on the right have zero cross-derivatives, they can be maximized separately. Let us 
set TT = (7ri,...,7rfc)', /3 = {l3[, . . . ,(3'^) , 0y == (0y|i, . . . , ^yn-) , Ox == (6'x|i, • • ■ , ^x|/c) and (f}x = 
(<j>x\i, • ■ • J '/'x ifc) • The maximum of equation (I14p with respect to tt, subject to the constraints on those 
parameters, is obtained by maximizing the augmented function 

n k / ^ 

i=l j=l \j=l 

where A is a Lagrangian multiplier. Setting the derivative of equation (jlSp with respect to ttj equal to 
zero and solving for ttj yields 

n 

Maximizing (J14p with respect to /3 (and eventually, to (py) is equivalent to independently maximizing 
each of the k expressions 

n 
i=l 

The maximization of (fT7|) is equivalent to the maximization problem of the generalized linear model 
(for the complete data), except that each observation {xi,yi) contributes to the log-likelihood for each 
component with a known weight r^' , which is obtained in the preceding E-step. Maximization of 
(IT71) . with respect to jSj (and eventually, to cfiyij if the adopted distri bution is not fixed in the e xpo- 



nential family) can be carried out numerically; details can be found in 



Wedel and DeSarbol (|l995[) and 



Wedel and Kamakural (J200ll pp. 120-124) since the generalized linear exponential CWM shares this 
part of complete-data likelihood with the finite mixture of generalized linear models discussed in these 
references. 

Finally, maximizing (J14p with respect to dx (and eventually (fix) is equivalent to independently 
maximizing each of the k expressions 



e^^ 



j) 



4=1 



\n[p{x^;ex\j,<Px\j)] ■ 



(18) 



The maximization of (|18|) is equivalent to the maximization problem of the exponential family (for the 
complete data), except that each observation Xi contributes to the log-likelihood for each component 
with a known weight t^'^ , which is obtained in the preceding E-step. Maximization of (|18p. with 
respect to Ox\j (and eventually 0xh) can be carried out, as before, numerically. 

In the framework of model-based clustering, the fitted model can be used to classify the n observations 
via the maximum a posteriori (MAP) classification induced by the final estimates (denoted, as usual, by 
hats) 



MAP(fy) = <^ 



1 if max^ {^j/i} occurs at component j 
otherwise 



(19) 







1, . . . , n and j = 1, . . . , fc. Note that the MAP classification is used in the analyses of Section [T] 



5 Computational issues 



Code for the EM algorithm described in Section|4]was written in the R computing environment (|R Development Core Tean 
20111 ). and it is available at |http: //www, economia.unict . it/punzo In its actual version, the algorithm 
considers the simple four distributions in Table [TJ they are able to cover four different data types, as 
highlighted in the second row of Table [TJ Hence, we obtain sixteen different bivariate distributions for 
(X, Y) . This highlights the flexibility of resulting approach and conflrms the potentiality of the gener- 
alized linear exponential CWM in the idea to consider all the existing exponential family distributions. 
The number d of free parameters in each of the sixteen models in displayed in Table [2j For the Gamma 
distribution note that, as motivated at the end of Section E?^ the "log" link is considered. 



EM initialization: Before running the EM algorithm, the choice of the starting values constitutes an 
imp ortant issue. The standard i nitialization consists in selecting a value for xjj^ ' . An alternative approach 
(see 



McLachlan and Peel 



20001 p. 54) consists in performing the first E-step by specifying, in equation 



Table 2: Number of free parameters for each of the sixteen generahzed hnear exponential CWMs. 



X 



Y 



Gaussian 



Gamma 



Poisson 



Binomial 



Gaussian 


6fc-l 


6k- 1 


5k- 1 


5k- 1 


Gamma 


6k ~1 


6k- 1 


5k- 1 


5k- 1 


Poisson 


5k -1 


5k- 1 


4fc- 1 


4fc- 1 


Binomial 


5k -1 


5k- 1 


4fc- 1 


4fc- 1 



(fOl). the values of r 



(1) 



Biernacki et al. 



2003 . 



\, ... ,71 and i — 1 , . . . , fc. Among the possible in itialization strategies (see 



Karlis and Xekalakill2003. and 



the R-package f lexmix (jLeischI 



2004 and 



Bagnato and Punzo . 

J 



2012l for details) - according to 



Griin and Leischll2008[ ) which allows to estimate finite mixtures 



of generalized linear models - a random initialization is repeated t times from different random positions 
and the solution maximizing the observed-data log-likelihood among these t runs is selected. In each run, 
the n vectors r^ — I r^^ , . . . , r^). j can be randomly generated in a "soft" way, that is by generating k 
positive values summing to one, or alternatively in a "hard" way by randomly drawn from a multinomial 
distribution with probabilities (l//c,...,l/fc) . 



Convergence criterion: The Aitken acceleration procedure (JAitkenl . Il926r) is used to estimate the 
asymptotic maximum of the log-likelihood at each iteration of the EM algorithm. Based on this estimate, 
a decision can be made regarding whether or not the algorithm has reached convergence; that is, whether 
or not the log-likelihood is sufficiently close to its estimated asymptotic value. The Aitken acceleration 
at iteration g + 1, g = 0, 1, . . ., is given by 



,(9+1) 



/(9+2) „ l(q+l) 

l(q+i} - l(q) '' 



where /'^''+^', Z^'^-'^', and l'^'^^ are the log-likelihood values from iteratio ns q -t- 2, q -I- 1, and q, respectively. 



Then, the asymptotic estimate of the log-likelihood at iteration q + 2 (JBohning et al. 



19941 ) is given by 



/(9+2) ^ ;(«+!) _^ 



1 



1 - a(9+i) 



/(9+2) _ ^(9+1) 



In the analyses in Section [71 we follow 
e, with e = 0.05. 



Subedi et al 



( 2013 ) and stop our algorithms when Zoo — Z^'^^-* < 



Standard errors of the estimates: Once the EM algorithm is run, the covariance matrix of xp is 
calculated using the inverted negative Hessian matrix. The Hessian matrix is computed by the function 



10 



hessian in the R-package numDeriv. In particular, hessian is evaluated on the observed-data log- 
likelihood in correspondence to the solution provided by the EM algorithm. 

6 Model selection and clustering performance 

The generalized linear exponential CWM, in addition to xp, is also characterized by the number of 
components k. So far, this quantity has been treated as a priori fixed. Nevertheless, for practical 
purposes, choosing a relevant model needs its choice. In model-based clustering, model selection criteria 
are commonly used with this end. Among them, we will adopt the Bayesian information criterion (BIC; 



Schwarz . 



mm 



BIC = 2l(xpj -dln(n). 
In order to evaluate the clu stering performance in ca ses in which the true classification is known. 



the adjusted Rand index (ARI; 



Hubert and Arable 



19851) . and the misclassification rate (proportion of 



misallocated observations), will be taken into account. We recall that the ARI has an expected value of 
and perfect classification would result in a value equal to 1. The ARI and the misclassification rate 
will be respectively computed according to the functions adjust edRandlndex and classError of the 
mclust package for R. 



7 Illustrative examples 

This section presents some applications of the proposed model on both artificial and real data. Compar- 
isons with finite mixtures of generalized linear models of equation ([5]) , and with exponential-exponential 
LCMs of equation (llOp . are also considered; the former are estimated via the stepFlexmix function of the 
R-package f lexmix, while a specific code has been implemented in R to fit the latter (always available 
at http: //www. economia.unict . it/punzo[ ). The number of repetitions, for the random initialization, 
will be fixed at i = 10 for all the considered models. 



7.1 Artificial data analysis 

This example is based on artificial data randomly generated from a Binomial-Poisson CWM (Binomial 
on Y\x,j and Poisson on X\j). The data consist of n = 300 bivariate observations with k ~ 2 groups, 
say Gi and G2, having size s ni = 712 — 150 . They are displ ayed in Figur e [1] by the CW-plot, a graphical 



representation proposed in 



Ingrassia et al 



(2013. see also 



Punzo . 



20121 ) to highlight some key aspects 



of data modeling through the CWM. Figure [T] (on the top) displays the bar plot of X, including the 
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Figure 1: CW-plot of the simulated data from a Binomial-Poisson CWM (n = 300 and k ^ 2). 
Coinciding points are marked by a single symbol only. 

overall empirical probability mass function as well as the empirical probability mass functions for Gi and 
G2 weighted by rij/n, j — 1,2, respectively; bars are color-coded, using a gray scale, according to the 
groups. The vertical segments, which are superimposed on each bar, allow to visualize the corresponding 
generating (true) probabilities. The scatter plot of the data, with the true Binomial regressions, is 
displayed at the bottom of Figure [TJ 

7.1.1 Clustering and fitting evaluation 

We aim at comparing the performance of model ([7]) with a finite mixture of Binomial regressions and 
with a Binomial-Poisson LCM (Binomial on Y\j and Poisson on X\j). The three models are directly 
fitted with fc = 2 and the corresponding scatter plots are shown in Figure[51 Here, it is straightforward to 
note as, although the true groups are well-separated (see the scatter plot in Figure [T]), only the Binomial- 
Poisson CWM produces good results. Table [3] also provides details about the obtained classifications 
which naturally confirm the considerations above. 

As concern the indirect applications of the models, apart from the mixture of Binomial regressions, the 
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(a) mixture of Binomial regressions (b) Binomial-Poisson LCM (fc = 2) (c) Binomial-Poisson CWM {k = 2) 

Figure 2: Estimated regression curves, and corresponding classification, for the NPreg data. 



indices 



models mixture of 

Binomial regressions 



Binomial-Poisson 
LCM 



Binomial-Poisson 
CWM 



ARI 
misclassification rate 



0.311 
0.220 



0.014 
0.437 



0.987 
0.003 



Table 3: Classification results, on the artificial data, using different mixture-based approaches (k = 2). 



remaining approaches are conceived to model the joint probability mass function of the data. Figure 3(a) 
shows the true probability mass function, while Figures |3(b)| and 3(c) display the fitted counterparts 
from the LCM and the CWM, respectively. Also in this case, only the Binomial-Poisson CWM is able 
to reproduce the underlying probability mass function. 






(a) True (b) Binomial-Poisson LCM (k = 2) (c) Binomial-Poisson CWM (k = 2) 

Figure 3: Bivariate probability mass function of the artificial data. 
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7.1.2 Classification evaluation 



Assume that one of the models introduced before has been fitted based on a sample A^ of size m. In 
this section we analyze the behavior of the fitted model in classifying new observations in some set 7", 
of size i, arising from the same statistical population. In the machine learning literature, M and T are 
called the training set and the test set, respectively. 

With this aim, a simple simulation study was performed by using the n = 300 artificial data of 
Section 17.11 The key aspect of the simulation study was to partition the 300 data into training and test 
sets so that n = rn + t = 300. In detail, for each value of m ranging from 150 to 270 (that is a proportion 
of the whole data ranging from 0.5 to 0.9), so that t ranged from 150 to 30 (that is a corresponding 
proportion of the whole data ranging from 0.5 to 0.1), five hundred vectors, of size m, were randomly 
generated from the set {1, . . . ,300} without replacement. Their elements indicated the observations to 
consider into the training set. In each of the five hundred replications, three models were fitted on the 
training set, directly with k — 2: the finite mixture of Binomial regressions, the Binomial-Poisson LCM, 
and the Binomial-Poisson CWM. Then, the m observations of the training set were classified by MAP 
while the t = 300 — m observations of the test set were classified, always by MAP, considering the fitted 
model. 

Figure |4] shows the average ARI values, computed across the 500 replications on both the training set 
(see Figure 4(a) ) and the test set (see Figure [4(b)[ ), for each considered value of to. Analogously, Figure[5] 
displays the average misclassification rates. Regardless from the considered index, and regardless from 
the set of observations on which the index is computed, the classification performance of the Binomial- 
Poisson CWM is always clearly the best, and it is very near to the optimality. This performance is also 
stable with respect to m. 



7.2 Real data analysis 



The real data of this example, also analyzed in 



Cicia et al. 



pOlOf l. are based on a sample oin = 224 fair- 



trade coffee consumers interviewed at stores. The variables of interest are: importance that respondents 
attribute to price in their fair-trade coffee purchase (Y*; measured on a 1-7 Likert scale where 1 indicated 
"completely unimportant" and 7 "extremely important") and number of fair-trade coffee packages over 
10 purchases {X*] tho se who did not con sume coffee were excluded from the survey leading to values 



ranging from 1 to 10). 



Cicia et al. 



(|2010r ) studied the probability mass function of Y* as function of a 



dichotomous version of X * assuming value 1 if X* < 4, and otherwise, by using the CUB model of 



D'Elia and Piccolo 



20051) . Our analysis aims at evaluating the validity of this dichotomization via the 



models described in Section [3] Thus, X* and Y* have been replaced by X = A* — 1 and Y = Y* - 1, 
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mixture of Binomial regressions 
Binomial-Poisson LCi^ 
Binomial-Poisson CWM 



mixture of Binomial regressions 
Binomial-Poisson LCM 
Binomial-Poisson CWM 



150 162 174 



= 10 222 234 246 258 270 



150 138 126 114 102 



54 42 30 



(a) training set Ai (b) test set T 

Figure 4: Average ARI values over 500 replications. 



mixture of Binomial regressions 
Binomial-Poisson LCtJ\ 
Binomial-Poisson CWM 



mixture of Binomial regres 
Binomial-Poisson LCM 
Binomial-Poisson CWM 



150 162 174 



198 210 222 234 246 258 270 



150 138 



114 102 90 78 



54 42 30 



(a) training set A4 (b) test set T 

Figure 5: Average misclassification rates over 500 replications. 

respectively, and a Binomial distribution is adopted for both X and Y^ with mx — 9 and my — 6 (see 
Table [I]). 

The CW-plot of the observed data, in absence of group membership information, is displayed in 
Figure[Sl Here, the size of each point in the scatter plot is proportional to the frequency of the pair {x. y) , 
for X — 0, 1, . . . , 9 and y = 0,l,...,6. This also gives information about the probability mass of the pairs. 
The models applied on these data are the Binomial-Binomial LCM and the Binomial-Binomial CWM. 
They are fitted with values of k ranging from 1 to 8. Results for the mixture of Binomial regressions are 
not reported here for two reasons: first, when A; > 2, the number of fitted mixture components is always 
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• • • • 



• 



• 



Figure 6: CW-plot of the real data {n = 224). In the scatter plot, the point size is proportional to the 
nunrbcr of occurrences. 

lower than k; second, when fc > 2, there is always a group having a very small number of pairs {x,y) 
inside. 

Figure [7] displays the BIG values for the two considered approaches. For both of them, the choice 
A; = 2 is the best one and a strong improvement, of the considered criterion, is obtained by moving 
from k ^ 1 to k = 2; this justifies the presence of groups in the data. The CW-plots of the best BIG 
models are given Figure |S1 The two approaches provide similar results, especially with reference to the 
regression model in the "black" group. Moreover, for bo th the models, the obtained classification is in 



agreement with the dichotomization of 



Cicia et al. 



( 201CI ). as also corroborated by the summary indices 



of Tabled A slight higher agreement is obtained for the Binomial-Binomial CWM. Finally recall that, 
due to Proposition[21 the Binomial-Binomial LCM can be seen as nested in the Binomial-Binomial CWM. 
This also means that the likelihoods, as well as the likelihood-based model selection criteria like the BIC, 
are comparable among models. Thus, based on the plot in Figure [71 we can select both the best number 
of mixture components for each model and the best model. In these terms, the best configuration is 
represented by the Binomial-Binomial LCM with k ^ 2: furthermore, apart from the case fc = 1, the 
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Figure 7: BIG values for the considered approaches 
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(a) Binomial-Poisson LCM (A: = 2) (b) Binomial-Poisson CWM {k = 2) 

Figure 8: Cw-plot of the fitted models on the real data. 

Binomial-Binomial LCM is always better than its CWM counterpart for each considered value of k. 
Thus, for these data, wc can conclude that: first, local indepe n dence can be assumed instead of local 



linear dependence; second, the dichotomization of 
mixtures with random covariates. 



Cicia et al 



pOld ) results justified by the adopted 
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models Binomial-Binomial Binomial-Binomial 

indices ^~--~---__^ j^^^ p^j^ 

ARI 0.749 0.812 

misclassification rate 0.067 0.049 

Table 4: Classification results, on the artificial data, using different mixture-based approaches {k — 2). 

8 Discussion and future work 



Ingrassia et al. 



( 2012a 



In the bivariate case, an extension of the generalized linear Gaussian CWM of 
was introduced to model covariates defined on various supports (e.g., finite, countable, and positive) as 
allowed by the exponential family distributions. The new model, named generalized linear exponential 
CWM, was justified on the ground of density estimation and model-based clustering and it was mainly 
motivated by the possibility to fit data of mixed-type. Theoretical arguments were also given showing the 
relationships with some well-known mixture-based approaches. The natural counterpart of the proposed 
model, in the frame of latent class modeling, was also given; it was named exponential-exponential LCM 
and was shown to be nested in the generalized linear Gaussian CWM. Parameter estimation was carried 
out within the EM algorithm framework, and the BIG was used for determining the number of groups. 
Moreover, given the relationship between the generalized linear Gaussian CWM and the exponential- 
exponential LCM, the BIG was also adopted to select among them; in other words, the BIG represents a 
way to evaluate if an exponential-exponential LCM, which assumes local independence, is good enough 
for data at hand. The proposed model was finally applied on real/artificial data where it shown a good 
performance. 

Future work will follow several avenues. To begin, extension of the generalized linear exponential 
CWM, to more than two variables, will be considered. This extension could initially concern only the X 
variable and then involve also the Y one. Furthermore, it could be interesting to evaluate, theoretically or 
by si mulation, the m o st con venient model selection criterion for the proposed model. This search, in line 



with 



Biernacki et al. 



(J200d . p. 724), could be separately conducted according to the type of application. 



direct or indirect, of the model. 

Proof of Proposition 1 

Given tt and ^y, if the component marginal densities of X do not depend from j, that la li 6 x\i ^ ■ ■ ■ = 
Ox\k — Ox and 4>x\i = ■ ■ ■ = 4>x\k = (f'x, then the posterior probabilities of component- membership for 
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the CWM in ([7]) can be written as 

TTjP [vi I Xi ■,f3j,(t)Y\j)p {xi ■,Ox,(t)x) 



P{Z,j =l\Tv,^Y^^x,(t)x) = ^, 



^ TTjp (yj \xi;f3j,(j)Y\j ) P (a;*; Ox,(l>x) 
TTjP (?/i I x.i ; /3j , 0y I J ) ^Xait^-^T^^ 

k 

T^jPivi \xi;f3j,(t)Y\j) 



k 
^■^jPivt \xi;(3j,(f>Y\j) 



i = 1, . . . , n and j = 1, . . . ,k, 



which coincide with the posterior probabihties associated to the model in (|S]). 

Proof of Proposition 2 

Given tt and ^x^ if th^ component generahzed hnear models does not depend from j, that is if /3i = 
• • ■ = /3j, = /3 and (/lyii = • ■ • = ^yij, = (/i)y, then the posterior probabilities of component- membership 
for model ([7]) can be written as 

p/7 -11 £ a A. \ TTjP {y, \xi; f3, (I>y)p{x,; 6>jc|j->x|j) 
F{Zij = 1\tt,^x^Pi9y) = — 

^ TTjp {yi \x^;|3,(t)Y)p {xi; Ox\j,(f>x\j) 
_ Trjp{xi;9x\j,(t)x\j) 



i = 1, . . . ,n and j = 1, . . . , fc, 



^7rjp(a;,;6'x|j,0x|j) 
which coincide with the posterior probabilities of the model in Q. 



Proof of Proposition 3 

The result comes from the one-to-one relation between Oyij and /3oj, j — 1, . . . , fc, which is due to the 
monotonicity of the link function. 
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